################################################
# Analysis for:
#
# Warner, Zach, Garrett N. Vande Kamp, and Soren Jordan.
# ``Conditional Relationships in Dynamic Models.'' 
#  Political Science Research & Methods
#
# General Supplemental Analysis (Section 5)
#
# Required files: WVKJ-MC-simulations.csv
#                 interp_urdf.R (for interpretation)
#
################################################

# Set your working directory here

library(dplyr)
library(ggplot2)
library(grid)
library(haven)
library(lmtest)
library(dynamac)
library(urca)
library(tseries)
library(polywog)
library(MASS)

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  if (is.null(layout)) {
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  if (numPlots==1) {
    print(plots[[1]])
  } else {
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    for (i in 1:numPlots) {
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

# ur.df output is annoying. interp_urdf.R helper function for interpretation
#  Author: Hank Roark. 
#  https://gist.github.com/hankroark/968fc28b767f1e43b5a33b151b771bf9
source("interp_urdf.R")

# These parameters are calso called in the sim analysis below the data creation.
# everything stationary with rho = 0.5
rho.x <- c(0.1, 0.5, 0.9) 
rho.z <- c(0.1, 0.5, 0.9) 
rho.y <- c(0.1, 0.5, 0.9) 
nobs <- c(16, 31, 51, 101, 501)  # add 1 because we'll drop the first value due to lagged values being NA

#set up the master dataframes for conditions to vary
exper <- expand.grid(nobs = nobs,
                     rho.x = rho.x,
                     rho.z = rho.z,
                     rho.y = rho.y)

# true values not varying
alpha.0 <- 0   # intercept
alpha.1 <- .5  # l.y
beta.0 <- -.3   # x
beta.1 <- .3   # l.x
beta.2 <- .2  # z
beta.3 <- -.7  # l.z
beta.4 <- .5  # x_z
beta.5 <- -.25  # l.x_z
beta.6 <- .1  # x_l.z
beta.7 <- -.4  # l.x_l.z


df <- read.csv("WVKJ-MC-simulations.csv")

df$var <- as.character(df$var)
df$model <- as.character(df$model)
df$dgp <- as.character(df$dgp)


### Reset the seed so we get the same results whether or not you've run the MC sims
set.seed(983405201)

# Add pretty plot labels
df$t <- df$nobs
df$t[df$t == 16] <- "T = 15"
df$t[df$t == 31] <- "T = 30"
df$t[df$t == 51] <- "T = 50"
df$t[df$t == 101] <- "T = 100"
df$t[df$t == 501] <- "T = 500"

df$t <- factor(df$t, levels = c("T = 15", "T = 30", "T = 50", "T = 100", "T = 500"))
table(df$nobs, df$t)

df$rho.y.fac <- recode_factor(df$rho.y, `0.1` = "rho[y] == 0.1", 
															`0.5` = "rho[y] == 0.5",
															`0.9` = "rho[y] == 0.9")

df$rho.x.fac <- recode_factor(df$rho.x, `0.1` = "rho[x] == 0.1", 
															`0.5` = "rho[x] == 0.5",
															`0.9` = "rho[x] == 0.9")

df$rho.z.fac <- recode_factor(df$rho.z, `0.1` = "rho[z] == 0.1", 
															`0.5` = "rho[z] == 0.5",
															`0.9` = "rho[z] == 0.9")

df$model.fac <- df$model
df$model.fac[df$model.fac == "full"] <- "Full interaction"
df$model.fac[df$model.fac == "restricted"] <- "Restricted"
df$model.fac[df$model.fac == "noint"] <- "No interaction"
df$model.fac <- factor(df$model.fac, levels = c("No interaction", "Restricted", "Full interaction"))


# Finally, define a df with reasonable ts
df.ts <- df[which(df$nobs %in% c(51, 101, 501)), ]






##############################################################################################################################
# Figure 1 in the main text averages over rho.x, rho.z, and t. let's see if that's a problem
##############################################################################################################################

### Problem 1: bias in ECR. We should see (1) worse estimates, (2) lower coverage rates
# alpha estimates
a1 <- df.ts[which(df.ts$var == "alpha.1"),]
a1$true <- ifelse(a1$dgp == "spurious", a1$rho.y, alpha.1)
a1$location <- as.numeric(factor(a1$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
a1$location[which(a1$model == "noint")] <- a1$location[which(a1$model == "noint")] - 0.1
a1$location[which(a1$model == "full")] <- a1$location[which(a1$model == "full")] + 0.1

# coverage
cov1 <- df.ts[which(df.ts$var == "cov.alpha.1"),]
cov1$location <- as.numeric(factor(cov1$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
cov1$location[which(cov1$model == "noint")] <- cov1$location[which(cov1$model == "noint")] - 0.1
cov1$location[which(cov1$model == "full")] <- cov1$location[which(cov1$model == "full")] + 0.1
cov1$true <- rep(.95, nrow(cov1))


## facet by all the rhos; separate for each T
# power: facet by t
fig1sm.1.t50 <- ggplot(a1[a1$nobs == 51,], aes(x = location, y = median, color = as.factor(rho.y))) +
  theme_bw() +
  facet_grid(vars(rho.x.fac), vars(rho.z.fac), labeller = label_parsed) + 
  geom_segment(aes(x = .75, xend = 1.25, y = rho.y, yend = rho.y), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") + 
  geom_segment(aes(x = 2.75, xend = 3.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = model.fac), size = I(1.8), alpha = I(0.5)) +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("L", "N", "R", "F")) +
  labs(title = "ECR estimates: T = 50", x = "", 
  	   color = expression(paste(rho[y])), 
       shape = "Model", y = expression(paste("Median ", hat(alpha)[1], " for each set of 500 replicates", sep = ""))) 

fig1sm.2.t50 <- ggplot(cov1[cov1$nobs == 51,], aes(x = location, y = mean,  color = as.factor(rho.y))) +
  theme_bw() +
  facet_grid(vars(rho.x.fac), vars(rho.z.fac), labeller = label_parsed) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = model.fac), size = I(1.8), alpha = I(0.5)) +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("L", "N", "R", "F")) +
  labs(title = "ECR coverage: T = 50", x = "Data-generating process: 'L'DV, 'N'o interaction,\n'R'estricted, or 'F'ull interaction",
  		color = expression(paste(rho[y])), 
       shape = "Model", y = expression(paste(hat(alpha)[1], " coverage probability for each set of 500 replicates", sep = ""))) 

######### Figure SM. 1 #########
pdf("sm-fig1-byrho-t50.pdf", width = 6, height = 8)
multiplot(fig1sm.1.t50, fig1sm.2.t50, cols = 1)
dev.off()


fig1sm.1.t100 <- ggplot(a1[a1$nobs == 101,], aes(x = location, y = median, color = as.factor(rho.y))) +
  theme_bw() +
  facet_grid(vars(rho.x.fac), vars(rho.z.fac), labeller = label_parsed) + 
  geom_segment(aes(x = .75, xend = 1.25, y = rho.y, yend = rho.y), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") + 
  geom_segment(aes(x = 2.75, xend = 3.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = model.fac), size = I(1.8), alpha = I(0.5)) +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1),  labels = c("L", "N", "R", "F")) +
  labs(title = "ECR estimates: T = 100", x = "",
		 color = expression(paste(rho[y])), 
       shape = "Model", y = expression(paste("Median ", hat(alpha)[1], " for each set of 500 replicates", sep = ""))) 

fig1sm.2.t100 <- ggplot(cov1[cov1$nobs == 101,], aes(x = location, y = mean,  color = as.factor(rho.y))) +
  theme_bw() +
  facet_grid(vars(rho.x.fac), vars(rho.z.fac), labeller = label_parsed) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = model.fac), size = I(1.8), alpha = I(0.5)) +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1),  labels = c("L", "N", "R", "F")) +
  labs(title = "ECR coverage: T = 100", x = "Data-generating process: 'L'DV, 'N'o interaction,\n'R'estricted, or 'F'ull interaction",
		color = expression(paste(rho[y])), 
       shape = "Model", y = expression(paste(hat(alpha)[1], " coverage probability for each set of 500 replicates", sep = ""))) 


######### Figure SM. 2 #########
pdf("sm-fig1-byrho-t100.pdf", width = 6, height = 8)
multiplot(fig1sm.1.t100, fig1sm.2.t100, cols = 1)
dev.off()

fig1sm.1.t500 <- ggplot(a1[a1$nobs == 501,], aes(x = location, y = median, color = as.factor(rho.y))) +
  theme_bw() +
  facet_grid(vars(rho.x.fac), vars(rho.z.fac), labeller = label_parsed) + 
  geom_segment(aes(x = .75, xend = 1.25, y = rho.y, yend = rho.y), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") + 
  geom_segment(aes(x = 2.75, xend = 3.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = model.fac), size = I(1.8), alpha = I(0.5)) +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1),  labels = c("L", "N", "R", "F")) +
  labs(title = "ECR estimates: T = 500", x = "",
		color = expression(paste(rho[y])), 
       shape = "Model", y = expression(paste("Median ", hat(alpha)[1], " for each set of 500 replicates", sep = ""))) 

fig1sm.2.t500 <- ggplot(cov1[cov1$nobs == 501,], aes(x = location, y = mean,  color = as.factor(rho.y))) +
  theme_bw() +
  facet_grid(vars(rho.x.fac), vars(rho.z.fac), labeller = label_parsed) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = model.fac), size = I(1.8), alpha = I(0.5)) +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("L", "N", "R", "F")) +
  labs(title = "ECR coverage: T = 500", x =  "Data-generating process: 'L'DV, 'N'o interaction,\n'R'estricted, or 'F'ull interaction",
		color = expression(paste(rho[y])), 
       shape = "Model", y = expression(paste(hat(alpha)[1], " coverage probability for each set of 500 replicates", sep = ""))) 


######### Figure SM. 3 #########
pdf("sm-fig1-byrho-t500.pdf", width = 6, height = 8)
multiplot(fig1sm.1.t500, fig1sm.2.t500, cols = 1)
dev.off()





##############################################################################################################################
# Figure 2 in the main text averages over rho.x, rho.z. let's see if that's a problem
##############################################################################################################################

rmse <- df.ts[which(df.ts$var == "RMSE"), ]
rmse$location <- as.numeric(factor(rmse$dgp, 
                                   levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
rmse$location[which(rmse$model == "noint")] <- 
  rmse$location[which(rmse$model == "noint")] - 0.1
rmse$location[which(rmse$model == "full")] <- 
  rmse$location[which(rmse$model == "full")] + 0.1


######### Figure SM. 4 #########
pdf("sm-fig2-byrho-t50.pdf", width = 6, height = 8)
ggplot(rmse[rmse$nobs == 51,], aes(x = location, y = mean,  color = as.factor(rho.y))) +
  theme_bw() +
  facet_grid(vars(rho.x.fac), vars(rho.z.fac), labeller = label_parsed) + 
  theme(panel.grid.minor = element_blank(), panel.border = element_rect(fill=NA, size=1, linetype = "solid")) +
  theme_bw() +
  geom_point(aes(shape = model.fac), size = I(1.8), alpha = I(0.5)) +
  theme(legend.position = "none") +
 # coord_cartesian(ylim = c(.98, 1.88)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("L", "N", "R", "F")) +
  labs(x =  "Data-generating process: 'L'DV, 'N'o interaction,\n'R'estricted, or 'F'ull interaction",
		y = "RMSE", color = expression(paste(rho[y])), 
       shape = "Model")
dev.off()


######### Figure SM. 5 #########
pdf("sm-fig2-byrho-t100.pdf", width = 6, height = 8)
ggplot(rmse[rmse$nobs == 101,], aes(x = location, y = mean,  color = as.factor(rho.y))) +
  theme_bw() +
  facet_grid(vars(rho.x.fac), vars(rho.z.fac), labeller = label_parsed) + 
  theme(panel.grid.minor = element_blank(), panel.border = element_rect(fill=NA, size=1, linetype = "solid")) +
  theme_bw() +
  geom_point(aes(shape = model.fac), size = I(1.8), alpha = I(0.5)) +
  theme(legend.position = "none") +
 # coord_cartesian(ylim = c(.98, 1.88)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("L", "N", "R", "F")) +
  labs(x =  "Data-generating process: 'L'DV, 'N'o interaction,\n'R'estricted, or 'F'ull interaction",
		y = "RMSE", color = expression(paste(rho[y])), 
       shape = "Model")
dev.off()


######### Figure SM. 6 #########
pdf("sm-fig2-byrho-t500.pdf", width = 6, height = 8)
ggplot(rmse[rmse$nobs == 501,], aes(x = location, y = mean,  color = as.factor(rho.y))) +
  theme_bw() +
  facet_grid(vars(rho.x.fac), vars(rho.z.fac), labeller = label_parsed) + 
  theme(panel.grid.minor = element_blank(), panel.border = element_rect(fill=NA, size=1, linetype = "solid")) +
  theme_bw() +
  geom_point(aes(shape = model.fac), size = I(1.8), alpha = I(0.5)) +
  theme(legend.position = "none") +
 # coord_cartesian(ylim = c(.98, 1.88)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("L", "N", "R", "F")) +
  labs(x =  "Data-generating process: 'L'DV, 'N'o interaction,\n'R'estricted, or 'F'ull interaction",
		y = "RMSE", color = expression(paste(rho[y])), 
       shape = "Model")
dev.off()




##############################################################################################################################
# Figure 1 in the main text only explores alpha.1. What about all of the beta estimates?
##############################################################################################################################
###################################################################################################
## Beta 0 (x)
b0 <- df.ts[which(df.ts$var == "beta.0"),]
b0$true <- ifelse(b0$dgp == "spurious", 0, beta.0)
b0$location <- as.numeric(factor(b0$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
b0$location[which(b0$model == "noint")] <- b0$location[which(b0$model == "noint")] - 0.1
b0$location[which(b0$model == "full")] <- b0$location[which(b0$model == "full")] + 0.1

# coverage
covb0 <- df.ts[which(df.ts$var == "cov.beta.0"),]
covb0$location <- as.numeric(factor(covb0$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
covb0$location[which(covb0$model == "noint")] <- covb0$location[which(covb0$model == "noint")] - 0.1
covb0$location[which(covb0$model == "full")] <- covb0$location[which(covb0$model == "full")] + 0.1
covb0$true <- rep(.95, nrow(covb0))


figsmb0.1 <- ggplot(b0, aes(x = location, y = median)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = beta.0, yend = beta.0), color = "gray30", linetype = "dashed") + 
  geom_segment(aes(x = 2.75, xend = 3.25, y = beta.0, yend = beta.0), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = beta.0, yend = beta.0), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
   theme(legend.position = "none") +
  # coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter estimates", x = "Data-generating process",
       y = expression(paste("Median ", hat(beta)[0], " for each set of 500 replicates", sep = ""))) 
       
figsmb0.2 <- ggplot(covb0, aes(x = location, y = mean)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter coverage", x = "Data-generating process",
       y = expression(paste(hat(beta)[0], " coverage probability for each set of 500 replicates", sep = "")))


######### Figure SM. 7 #########
pdf("sm-b0.pdf", width = 6, height = 8)
multiplot(figsmb0.1, figsmb0.2, cols = 1)
dev.off()



###################################################################################################
## Beta 1 (x[t-1])
b1 <- df.ts[which(df.ts$var == "beta.1"),]
b1$true <- ifelse(b1$dgp == "spurious", 0, beta.1)
b1$location <- as.numeric(factor(b1$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
b1$location[which(b1$model == "noint")] <- b1$location[which(b1$model == "noint")] - 0.1
b1$location[which(b1$model == "full")] <- b1$location[which(b1$model == "full")] + 0.1

# coverage
covb1 <- df.ts[which(df.ts$var == "cov.beta.1"),]
covb1$location <- as.numeric(factor(covb1$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
covb1$location[which(covb1$model == "noint")] <- covb1$location[which(covb1$model == "noint")] - 0.1
covb1$location[which(covb1$model == "full")] <- covb1$location[which(covb1$model == "full")] + 0.1
covb1$true <- rep(.95, nrow(covb1))


figsmb1.1 <- ggplot(b1, aes(x = location, y = median)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = beta.1, yend = beta.1), color = "gray30", linetype = "dashed") + 
  geom_segment(aes(x = 2.75, xend = 3.25, y = beta.1, yend = beta.1), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = beta.1, yend = beta.1), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
   theme(legend.position = "none") +
  # coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter estimates", x = "Data-generating process",
       y = expression(paste("Median ", hat(beta)[1], " for each set of 500 replicates", sep = ""))) 
       
figsmb1.2 <- ggplot(covb1, aes(x = location, y = mean)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter coverage", x = "Data-generating process",
       y = expression(paste(hat(beta)[1], " coverage probability for each set of 500 replicates", sep = "")))


######### Figure SM. 8 #########
pdf("sm-b1.pdf", width = 6, height = 8)
multiplot(figsmb1.1, figsmb1.2, cols = 1)
dev.off()



###################################################################################################
## Beta 2 (z)
b2 <- df.ts[which(df.ts$var == "beta.2"),]
b2$true <- ifelse(b2$dgp == "spurious", 0, beta.2)
b2$location <- as.numeric(factor(b2$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
b2$location[which(b2$model == "noint")] <- b2$location[which(b2$model == "noint")] - 0.1
b2$location[which(b2$model == "full")] <- b2$location[which(b2$model == "full")] + 0.1

# coverage
covb2 <- df.ts[which(df.ts$var == "cov.beta.2"),]
covb2$location <- as.numeric(factor(covb2$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
covb2$location[which(covb2$model == "noint")] <- covb2$location[which(covb2$model == "noint")] - 0.1
covb2$location[which(covb2$model == "full")] <- covb2$location[which(covb2$model == "full")] + 0.1
covb2$true <- rep(.95, nrow(covb2))


figsmb2.1 <- ggplot(b2, aes(x = location, y = median)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = beta.2, yend = beta.2), color = "gray30", linetype = "dashed") + 
  geom_segment(aes(x = 2.75, xend = 3.25, y = beta.2, yend = beta.2), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = beta.2, yend = beta.2), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
   theme(legend.position = "none") +
  # coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter estimates", x = "Data-generating process",
       y = expression(paste("Median ", hat(beta)[2], " for each set of 500 replicates", sep = ""))) 
       
figsmb2.2 <- ggplot(covb2, aes(x = location, y = mean)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter coverage", x = "Data-generating process",
       y = expression(paste(hat(beta)[2], " coverage probability for each set of 500 replicates", sep = "")))


######### Figure SM. 9 #########
pdf("sm-b2.pdf", width = 6, height = 8)
multiplot(figsmb2.1, figsmb2.2, cols = 1)
dev.off()



###################################################################################################
## Beta 3 (z[t-1])
b3 <- df.ts[which(df.ts$var == "beta.3"),]
b3$true <- ifelse(b3$dgp == "spurious", 0, beta.3)
b3$location <- as.numeric(factor(b3$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
b3$location[which(b3$model == "noint")] <- b3$location[which(b3$model == "noint")] - 0.1
b3$location[which(b3$model == "full")] <- b3$location[which(b3$model == "full")] + 0.1

# coverage
covb3 <- df.ts[which(df.ts$var == "cov.beta.3"),]
covb3$location <- as.numeric(factor(covb3$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
covb3$location[which(covb3$model == "noint")] <- covb3$location[which(covb3$model == "noint")] - 0.1
covb3$location[which(covb3$model == "full")] <- covb3$location[which(covb3$model == "full")] + 0.1
covb3$true <- rep(.95, nrow(covb3))


figsmb3.1 <- ggplot(b3, aes(x = location, y = median)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = beta.3, yend = beta.3), color = "gray30", linetype = "dashed") + 
  geom_segment(aes(x = 2.75, xend = 3.25, y = beta.3, yend = beta.3), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = beta.3, yend = beta.3), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
   theme(legend.position = "none") +
  # coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter estimates", x = "Data-generating process",
       y = expression(paste("Median ", hat(beta)[3], " for each set of 500 replicates", sep = ""))) 
       
figsmb3.2 <- ggplot(covb3, aes(x = location, y = mean)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter coverage", x = "Data-generating process",
       y = expression(paste(hat(beta)[3], " coverage probability for each set of 500 replicates", sep = "")))


######### Figure SM. 10 #########
pdf("sm-b3.pdf", width = 6, height = 8)
multiplot(figsmb3.1, figsmb3.2, cols = 1)
dev.off()



###################################################################################################
## Beta 4 (x*z)
b4 <- df.ts[which(df.ts$var == "beta.4"),]
b4$true <- ifelse(b4$dgp %in% c("spurious", "no_interaction"), 0, beta.4)
b4$location <- as.numeric(factor(b4$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
b4$location[which(b4$model == "noint")] <- b4$location[which(b4$model == "noint")] - 0.1
b4$location[which(b4$model == "full")] <- b4$location[which(b4$model == "full")] + 0.1

# coverage
covb4 <- df.ts[which(df.ts$var == "cov.beta.4"),]
covb4$location <- as.numeric(factor(covb4$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
covb4$location[which(covb4$model == "noint")] <- covb4$location[which(covb4$model == "noint")] - 0.1
covb4$location[which(covb4$model == "full")] <- covb4$location[which(covb4$model == "full")] + 0.1
covb4$true <- rep(.95, nrow(covb4))


figsmb4.1 <- ggplot(b4, aes(x = location, y = median)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") + 
  geom_segment(aes(x = 2.75, xend = 3.25, y = beta.4, yend = beta.4), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = beta.4, yend = beta.4), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
   theme(legend.position = "none") +
  # coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter estimates", x = "Data-generating process",
       y = expression(paste("Median ", hat(beta)[4], " for each set of 500 replicates", sep = ""))) 
       
figsmb4.2 <- ggplot(covb4, aes(x = location, y = mean)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter coverage", x = "Data-generating process",
       y = expression(paste(hat(beta)[4], " coverage probability for each set of 500 replicates", sep = "")))


######### Figure SM. 11 #########
pdf("sm-b4.pdf", width = 6, height = 8)
multiplot(figsmb4.1, figsmb4.2, cols = 1)
dev.off()



###################################################################################################
## Beta 5 (x[t-1]*z)
b5 <- df.ts[which(df.ts$var == "beta.5"),]
b5$true <- ifelse(b5$dgp == "full_interaction", beta.5, 0)
b5$location <- as.numeric(factor(b5$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
b5$location[which(b5$model == "noint")] <- b5$location[which(b5$model == "noint")] - 0.1
b5$location[which(b5$model == "full")] <- b5$location[which(b5$model == "full")] + 0.1

# coverage
covb5 <- df.ts[which(df.ts$var == "cov.beta.5"),]
covb5$location <- as.numeric(factor(covb5$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
covb5$location[which(covb5$model == "noint")] <- covb5$location[which(covb5$model == "noint")] - 0.1
covb5$location[which(covb5$model == "full")] <- covb5$location[which(covb5$model == "full")] + 0.1
covb5$true <- rep(.95, nrow(covb5))


figsmb5.1 <- ggplot(b5, aes(x = location, y = median)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") + 
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = beta.5, yend = beta.5), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
   theme(legend.position = "none") +
  # coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter estimates", x = "Data-generating process",
       y = expression(paste("Median ", hat(beta)[5], " for each set of 500 replicates", sep = ""))) 
       
figsmb5.2 <- ggplot(covb5, aes(x = location, y = mean)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter coverage", x = "Data-generating process",
       y = expression(paste(hat(beta)[5], " coverage probability for each set of 500 replicates", sep = "")))


######### Figure SM. 12 #########
pdf("sm-b5.pdf", width = 6, height = 8)
multiplot(figsmb5.1, figsmb5.2, cols = 1)
dev.off()


###################################################################################################
## Beta 6 (x*z[t-1])
b6 <- df.ts[which(df.ts$var == "beta.6"),]
b6$true <- ifelse(b6$dgp == "full_interaction", beta.6, 0)
b6$location <- as.numeric(factor(b6$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
b6$location[which(b6$model == "noint")] <- b6$location[which(b6$model == "noint")] - 0.1
b6$location[which(b6$model == "full")] <- b6$location[which(b6$model == "full")] + 0.1

# coverage
covb6 <- df.ts[which(df.ts$var == "cov.beta.6"),]
covb6$location <- as.numeric(factor(covb6$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
covb6$location[which(covb6$model == "noint")] <- covb6$location[which(covb6$model == "noint")] - 0.1
covb6$location[which(covb6$model == "full")] <- covb6$location[which(covb6$model == "full")] + 0.1
covb6$true <- rep(.95, nrow(covb6))


figsmb6.1 <- ggplot(b6, aes(x = location, y = median)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") + 
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = beta.6, yend = beta.6), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
   theme(legend.position = "none") +
  # coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter estimates", x = "Data-generating process",
       y = expression(paste("Median ", hat(beta)[6], " for each set of 500 replicates", sep = ""))) 
       
figsmb6.2 <- ggplot(covb6, aes(x = location, y = mean)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter coverage", x = "Data-generating process",
       y = expression(paste(hat(beta)[6], " coverage probability for each set of 500 replicates", sep = "")))


######### Figure SM. 13 #########
pdf("sm-b6.pdf", width = 6, height = 8)
multiplot(figsmb6.1, figsmb6.2, cols = 1)
dev.off()


###################################################################################################
## Beta 7 (x[t-1]*z[t-1])
b7 <- df.ts[which(df.ts$var == "beta.7"),]
b7$true <- ifelse(b7$dgp == "full_interaction", beta.7, 0)
b7$location <- as.numeric(factor(b7$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
b7$location[which(b7$model == "noint")] <- b7$location[which(b7$model == "noint")] - 0.1
b7$location[which(b7$model == "full")] <- b7$location[which(b7$model == "full")] + 0.1

# coverage
covb7 <- df.ts[which(df.ts$var == "cov.beta.7"),]
covb7$location <- as.numeric(factor(covb7$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
covb7$location[which(covb7$model == "noint")] <- covb7$location[which(covb7$model == "noint")] - 0.1
covb7$location[which(covb7$model == "full")] <- covb7$location[which(covb7$model == "full")] + 0.1
covb7$true <- rep(.95, nrow(covb7))


figsmb7.1 <- ggplot(b7, aes(x = location, y = median)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") + 
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0, yend = 0), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = beta.7, yend = beta.7), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
   theme(legend.position = "none") +
  # coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter estimates", x = "Data-generating process",
       y = expression(paste("Median ", hat(beta)[7], " for each set of 500 replicates", sep = ""))) 
       
figsmb7.2 <- ggplot(covb7, aes(x = location, y = mean)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "Parameter coverage", x = "Data-generating process",
       y = expression(paste(hat(beta)[7], " coverage probability for each set of 500 replicates", sep = "")))


######### Figure SM. 14 #########
pdf("sm-b7.pdf", width = 6, height = 8)
multiplot(figsmb7.1, figsmb7.2, cols = 1)
dev.off()





##############################################################################################################################
# Questions of power
##############################################################################################################################
# How often do we achieve significance?
stat.signs <- df.ts[which(df.ts$var %in% c("ss.beta.0", "ss.beta.1", "ss.beta.2", "ss.beta.3",
						"ss.beta.4", "ss.beta.5", "ss.beta.6", "ss.beta.7")),]
stat.signs.spurious <- stat.signs[which(stat.signs$dgp == "spurious"),]
stat.signs.spurious.full <- stat.signs.spurious[which(stat.signs.spurious$model == "full"),]
stat.signs.full <- stat.signs[which(stat.signs$model == "full"),]
stat.signs.full$dgp.fac <- NA
stat.signs.full$dgp.fac[stat.signs.full$dgp == "spurious"] <- "LDV"
stat.signs.full$dgp.fac[stat.signs.full$dgp == "no_interaction"] <- "No interaction"
stat.signs.full$dgp.fac[stat.signs.full$dgp == "restricted"] <- "Restricted"
stat.signs.full$dgp.fac[stat.signs.full$dgp == "full_interaction"] <- "Full interaction"

stat.signs.full$dgp.fac <- factor(stat.signs.full$dgp.fac, levels = c("LDV", "No interaction", "Restricted", "Full interaction"))


######### Figure SM. 15 #########
pdf("sm-power.pdf", width = 6, height = 5)
ggplot(stat.signs.full, aes(x = var, y = mean, color = dgp.fac)) +
  theme(panel.grid.minor = element_blank(), panel.border = element_rect(fill=NA, size=1, linetype = "solid")) +
  theme_bw() +
  facet_grid(vars(t)) + #, vars(var)) + 
  geom_boxplot() +
 # theme(legend.position = "none") +
 # coord_cartesian(ylim = c(.98, 1.88)) +
 scale_x_discrete(labels = expression(atop(paste(beta[0]), paste("(", x[t], ")")), 
 										atop(paste(beta[1]), paste("(", x[t-1], ")")),
 										atop(paste(beta[2]), paste("(", z[t], ")")),
 										atop(paste(beta[3]), paste("(", z[t-1], ")")),
 										atop(paste(beta[4]), paste("(", x[t], z[t], ")")),
 										atop(paste(beta[5]), paste("(", x[t-1], z[t], ")")),
 										atop(paste(beta[6]), paste("(", x[t], z[t-1], ")")),
 										atop(paste(beta[7]), paste("(", x[t-1], z[t-1], ")")))) +
  labs(x = "Coefficient",
       y = "Mean proportion of statistically significant estimates\nfor each set of 500 replicates",
       color = "DGP") #+
 # theme(text = element_text(family = "minpro"))
dev.off()


## ANother angle: is this model inappropriately significant for the spiurious dgp?
interaction.betas <- df.ts[which(df.ts$var %in% c("ss.beta.4", "ss.beta.5", "ss.beta.6", "ss.beta.7") & 
					df.ts$dgp == "spurious" &
					df.ts$model == "full"), ]

pdf("fig-statsign.pdf", width = 6, height = 3.5)
ggplot(interaction.betas, aes(x = var, y = mean)) +
  theme(panel.grid.minor = element_blank(), panel.border = element_rect(fill = NA, size = 1, linetype = "solid")) +
  theme_bw() +
  facet_grid(vars(t)) +
  geom_boxplot() +
  #theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  scale_x_discrete(labels = expression(paste(beta[4], " (", x[t], z[t], ")"),
  										paste(beta[5], " (", x[t-1], z[t], ")"),
  										paste(beta[6], " (", x[t], z[t-1], ")"),
  										paste(beta[7], " (", x[t-1], z[t-1], ")"))) +
  labs(title = NULL, x = "Coefficient", y = "Mean proportion of statistically significant\nestimates for each set of 500 replicates")# +
#  theme(text = element_text(family = "minpro"))
dev.off()

# No. (Unreported in SM)

